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ABSTRACT 

A variety of modeling techniques is used with surface photometry from the litera- 
ture and recently acquired high-accuracy stellar kinematic data to constrain the three- 
dimensional mass distribution in the luminous cuspy elliptical galaxy NGC 1700. First, 
we model the radial velocity field and photometry, and, using a Bayesian technique, es- 
timate the triaxiality T and short-to-long axis ratio c in five concentric annuli between 
approximately 1 and 3 effective radii. The results are completely consistent with T 
being constant inside about 2.5r e (36"; 6.7/i~ 1 kpc). Adding an assumption of constant 
T as prior information gives an upper limit of T < 0.16 (95% confidence); this relaxes 
to T < 0.22 if it is also assumed that there is perfect alignment between the angular 
momentum and the galaxy's intrinsic short axis. Near axisymmetry permits us then 
to use axisymmetric models to constrain the radial mass profile. Using the Jeans (mo- 
ment) equations, we demonstrate that 2-integral, constant-M/L models cannot fit the 
data; but a 2-integral model in which the cumulative enclosed M/L increases by a fac- 
tor of ~ 2 from the center out to 12/i -1 kpc can. Three-integral models constructed by 
quadratic programming show that, in fact, no constant-M/L model is consistent with 
the kinematics. Anisotropic 3-integral models with variable M/L, while not uniquely 
establishing a minimum acceptable halo mass, imply, as do the moment models, a cu- 
mulative M/Lb ~ 10/i at 12/i _1 kpc. We conclude that NGC 1700 represents the best 
stellar dynamical evidence to date for dark matter in elliptical galaxies. 

Subject headings: galaxies: elliptical and lenticular, cD — galaxies: individual (NGC 
1700) — galaxies: kinematics and dynamics — galaxies: structure 



Observations reported in this paper were obtained at the Multiple Mirror Telescope Observatory, a joint facility 
of the University of Arizona and the Smithsonian Institution. 
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1. Introduction 

Evidence is accumulating for a dichotomy between two classes of elliptical galaxies: high- 
luminosity, slowly rotating systems with shallow central brightness cusps and a tendency for boxy 
isophotes, and lower luminosity, rotationally supported systems with steeper central cusps and 
a tendency for diskiness (Lauer et al. 1995, Kormendy & Bender 1996, Faber et al. 1997, and 
references therein). Photometric shape differences between high- and low-luminosity ellipticals 
(Tremblay & Merritt 1996) are roughly consistent with the general expectation that rapidly rotating 
systems should have oblate axisymmetry while slow rotators are more likely triaxial. There are 
exciting suggestions that central cusps may drive evolution towards axisymmetry through orbital 
stochasticity (Merritt & Fridman 1996, Merritt & Valluri 1996). But where the causal arrows 
actually point among these various properties is still far from clear. 

There is consequently a growing need for intrinsic shape measurements for a subsample of well 
studied ellipticals. This paper continues our efforts to obtain such measurements. Previously we 
have derived approximate shapes for the "standard" elliptical NGC 3379 (Statler 1994c; hereafter 
S94c) and NGC 4589 (Statler 1994b; S94b) using photometry and kinematic data from the liter- 
ature. In this paper we take up the case of the high surface brightness elliptical NGC 1700. The 
availability of high quality multi-position-angle kinematics allows an intrinsic shape determination 
of unprecedented precision, as well as interesting constraints on the presence of dark matter in the 
outer parts of the system. 

NGC 1700 is an elliptical galaxy of fairly average luminosity (Mb = —22.3 for Hq = 50kms _1 
Mpc^ 1 ), but comparatively high central dispersion (~ 230kms _1 ) and small (14") effective radius. 
This puts it roughly 2a off the Fundamental Plane in the sense of having uncharacteristically low 
M/L for objects of this mass. NGC 1700 is notable for being the most luminous of the steep-cusped 
"power law" ellipticals in the sample observed by HST (Lauer et al. 1995; most recent compilation 
in Faber et al. 1997). The high surface brightness at r e makes it particularly attractive for a study of 
intrinsic shape and M/L gradient, since long-slit spectra with good signal-to-noise can be obtained 
out to several optical scale lengths. 

Statler, Smecker-Hane, & Cecil (1996, hereafter SSC) present kinematic profiles for NGC 1700 
out to ~ 4r e on four position angles. This is the first data set to meet the requirements for accurate 
shape fitting spelled out in S94b; accuracy in the mean velocity field (VF) is better than lOkms -1 
inside 2r e and better than 15% of the maximum inside 3r e . In addition to the usual profiles of v, 
a, /13, and SSC derive profiles of the velocity field asymmetry parameters (Vi, V^Vjj) (Statler 
Sz Fry 1994), and show all three to be consistent with zero over a factor ~ 10 in radius. From the 
VF's high degree of symmetry, SSC surmise that the triaxiality of NGC 1700 is probably small; 
i.e., that the shape is close to oblate. They also argue qualitatively that the shallow slope of the 
projected RMS velocity profile is consistent with a standard dark halo. Our objective in this paper 
is to make both of these statements quantitative. 

We proceed as follows: in Section 2 we apply the dynamical shape fitting method of S94b 
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to SSC's kinematic data and photometry from Franx et al. (1989). We find that NGC 1700 is 
very nearly oblate. This result lets us exploit existing tools for modeling axisymmetric systems to 
assess the mass-to-light ratio. In § 3 we demonstrate that two-integral models with constant M/L 
are strongly inconsistent with the data, while a model with M/L increasing outward by a factor 
of ~ 2 over the range of the observations is viable. We then construct anisotropic three- integral 
models, and show that even anisotropy cannot completely obviate the need for an M/L gradient. 
We summarize and discuss the broader implications of these results in § 4. 



The shape fitting method is described in detail in S94b and S94c. The underlying dynamical 
models (Statler 1994a, hereafter S94a) rest on the following assumptions: (1) mean rotation in the 
galaxy arises from internal streaming in a potential with negligible figure rotation; (2) the mean- 
motion streamlines of short-axis and long-axis tube orbits can be represented by coordinate lines 
in a confocal coordinate system;0 (3) the luminosity density pi is stratified on similar ellipsoids, 
i.e., pL(r,6,(f>) = p~L(r)p* L (9,(f)); and (4) the velocity field obeys a "similar flow" ansatz outside the 
tangent point for a given line of sight, v(r, 6, (p) = v(r)v*(9, <p). 

Assumptions (3) and (4) are used only for projecting the models. As long as and 
P~L(r)v(r) decrease faster than r -2 , nearly all of the contribution to the projection integrals will 
come from radii near the tangent point, and the result will be insensitive to the structure at larger 
r. We refer to this approximation as quasi-local projection. The requirement on pL{ r )v{r) generally 
limits the validity of the models to regions where the rotation curve is not steeply rising, which for 
NGC 1700 means r ^ r e . 

A single model is described by the parameters (T, cl, ^, d), where T is the triaxiality of the total 
mass distribution, cl is the short-to-long axis ratio of the luminosity distribution, f2 = (9e,4>e) is 
the orientation of the galaxy relative to the line of sight, and the vector d represents the remaining 
dynamical parameters which are defined in S94b. The projected model predicts the ellipticity and 
the radial velocities on each sampled position angle, from which we calculate the likelihood of the 
observations. Repeating for ~ 10 7 models covering the full parameter space, we determine the 
multidimensional likelihood L(T, cl, d). The Bayesian estimate of the galaxy's shape is then the 
probability distribution P(T,cl), obtained by multiplying the likelihood by a model for the parent 
distribution from which Nature drew the galaxy and integrating over all other parameters: 



2. Intrinsic Shape 



2.1. 



Method 




(1) 



2 Anderson & Statler (1998) demonstrate the validity of this assumption for realistic potentials and show that the 
triaxiality parameter recovered from the streamlines closely agrees with the triaxiality of the mass distribution. 
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In Bayesian terminology, the integrand is called the posterior density and F* the prior. The 
marginal posterior density P is the integral of the posterior over the nuisance parameters. We use 
P to stand for any marginal density; for instance, the joint distribution in shape and orientation is 
given by 

p(T, CL ,n) = J ddF;(T,c L ,n,d)L(T lCL ,n,d). (2) 

We can also factor the parent distribution according to 

F;(T,c L ,Q,d) = ±-F p (T,c L )F$(d). (3) 

The first factor describes a population with random orientations^ and the dependence on d can 
be factored out without loss of generality because the d themselves can implicitly be functions of 
T and cl- 

In modeling NGC 3379 and NGC 4589 it was necessary to average the kinematics over a 
single large radial bin on each PA, yielding only a single average shape for each system. For 
NGC 1700 we can extend the method because we have excellent data over nearly a factor of 4 in 
radius in the region consistent with the assumptions of the models. SSC's radial binning gives 5 
independent measurements along each PA, with bin centers at 12'/6, 15'/6, 21'/0, 3C/6, and 46'/8. 
We first model each annulus independently in the quasi-local approximation. Then, assuming 
that the intrinsic principal axes remain aligned through the galaxy, we can combine the posterior 
probabilities requiring that the five annuli match their local kinematics and ellipticities in the same 
orientation. If we further assume that the triaxialities of the mass and light distributions are equal, 
we can add the constraint that the triaxiality profile must reproduce the observed isophotal twist. 
This procedure is described more fully in Appendix A. 

We argued previously (SSC) that the parallelogram-shaped isophotal distortions seen beyond 
R 60" arise from an incompletely phase-mixed, differentially precessing ring or disk, presumably 
acquired by accretion. If so, then the photometric and kinematic twists that start near R ~ 40" 
may be the result of intrinsic twisting of the isodensity surfaces rather than a shape gradient in 
an intrinsically aligned system. In fact, we strongly suspect that this is the case; nonetheless, in 
this paper we proceed under the assumption that the isodensity surfaces are aligned and derive 
constraints on the shape gradient in that context. If the outermost annulus is actually misaligned, 
our results for the inner four annuli still hold and are virtually unaffected, as we have verified by 
recomputing the posterior densities excluding the last annulus. 



2.2. Data 



Actually, the population of ellipticals may have a preferred orientation if, for instance, a significant number of 
E's are misclassified SO's seen far from edge-on. 
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2.2.1. Photometry 

Ellipticity and position angle profiles are taken from Franx et al. (1989) and averaged over 
the approximate ranges of SSC's radial bins. At smaller radii the intervals over which the profiles 
are averaged overlap slightly in order to better reflect the ellipticity gradient inside 25"; this has a 
negligible effect on the adopted profile but (intentionally) increases the error bars where the gradient 
is steeper. The results are given in columns 2 and 3 of Table |l|. Data are weighted according to 
the inverse square of Franx et a/.'s error bars. The tabulated uncertainty is the larger of the mean 
observational error per data point and the variance over the bin. 

The logarithmic slope of the surface brightness profile wavers non-monotonically between about 
— 1.7 and —2.3 for 10" <^ r <J 55". For simplicity in projection we adopt a power-law radial profile 
for the luminosity density, p ~ r~ k , compute models for k = (^-,3, ^p), and average the results. It 
turns out that the results are insensitive to the choice of k. 



2.2.2. Kinematics 

Mean velocities, corrected for the non-Gaussian shape of the LOSVD, are taken directly from 
columns 10 and 11 of SSC's Table 2 (as plotted in their Figure 10), folded about the origin, and 
averaged. Results are given in columns 4 - 7 of Table [l] for PAs (near minor axis), 225, 270, and 
315. 

The kinematic sampling is at a different position angle with respect to the photometric axes in 
each annulus because the major axis PA is not strictly constant in radius. However, the PA twist 
for 10" < r < 35" is so small that we assume the sampling in the inner four annuli is the same 
for the purpose of calculating the projected model velocities; this saves a factor of 3 in computing 
time. 

Over the region being modeled, the observed mean rotation amplitude is constant or slowly 
increasing with radius, though on some PAs the rotation curve does fall at large R. Again for 



Table 1. Photometric and Kinematic Data Used in Shape Fitting 



n 


c 




PA 


I'd 


V225 


«270 


«315 


12.6 


0.262 ± 


.006 


89.5 ± 1.0 


-5.6 ±13.1 


57.7 ± 5.5 


94.8 ± 6.0 


54.9 ± 8.4 


15.6 


0.270 ± 


.010 


89.0 ± 1.0 


-4.6 ± 6.5 


73.5 ± 5.6 


104.0 ± 6.1 


71.8 ± 10.4 


21.0 


0.283 ± 


.004 


89.0 ± 1.0 


12.5 ± 9.7 


79.2 ± 6.9 


113.8 ± 6.8 


71.3 ± 7.1 


30.6 


0.292 ± 


.003 


89.8 ±1.0 


-9.7 ± 11.9 


77.9 ± 10.3 


113.0 ± 9.8 


80.3 ± 10.8 


46.8 


0.296 ± 


.010 


95.8 ±3.2 


-10.6 ±30.3 


53.1 ±26.9 


133.8 ±15.0 


94.6 ± 17.4 
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simplicity in projection, we take a power-law form for the velocity scaling law v(r) ~ r~ l and 
compute models for I = (0, ±^). As for the density index k, the results are not sensitive to the 
value of I, and we simply average the derived probability densities. 

2.2.3. Forbidden Velocities on PA 225? 

SSC reported a striking reversal in the rotation profile along PA 225, at r ^ 40" to the 
northeast, and argued that this feature, most prominent in the outermost data point, is real. Two 
of us (TSS and TS-H) subsequently obtained additional spectra along this PA with the MMT Red 
Channel on 1995 February 3-4 UT. The instrument was rotated 180° relative to the earlier setup, 
and the galaxy was displaced ~ 60" to the southwest, putting the region of interest near the center 
of the slit. Reductions were performed as described by SSC. 

The resulting rotation curve shows no velocity reversal. We now believe the earlier result 
to have been spurious, caused by a distortion of the cross-correlation peak by some unnoticed or 
inadequately corrected systematic effect (possibly scattered light) in the weakly exposed outer part 
of the image. 

The data were reduced, unfortunately, after the present models had all been computed, but 
the results for intrinsic shape are not compromised. The radius at which the velocity had appeared 
to change sign is outside the region we are modeling. The revised rotation profile alters only our 
outermost annulus, through one data point which is averaged with its counterpart on the opposite 
side of the galaxy. Were we to include the revision, the value of i>225 in the last row of Table [l] 
would change from 53.1±26.9kms _1 to approximately 84±21kms~ 1 . This is a change of just over 
la in one datum (with the second-largest error bar) out of 20 in the velocity field. We have rerun 
models over part of the parameter space with this change included and find only small differences; 
we will point out below where the effect is noticeable. 

2.3. Results 

Models are computed over essentially the same grid of dynamical parameters used previously 
(S94b, S94c). At each of the 9 pairs of (k,l) values there are 8 choices of the "contrast function" 
(giving the relative streaming amplitude in short-axis and long-axis tubes), 4 choices of the function 
v*(t) (giving the angular dependence of the mean velocity away from the symmetry planes), and 
two treatments of radial mean motions, with streamlines assumed to lie on either spherical or 
ellipsoidal shells. The roles of the contrast and v*(t) functions are summarized in Appendix B. We 
do not calculate models in which the mass and luminosity distributions have different triaxialities 
since our earlier tests showed that the triaxiality of the latter has little bearing on that of the 
former. When we include the alignment and isophotal twist constraints we implicitly assume the 
two triaxialities are equal; this is justifiable a posteriori since nearly axisymmetric potentials turn 
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out strongly preferred. By Monte Carlo sampling over (T, c, the likelihood L(T, c, f2, d) for 
each annulus is computed for each set of dynamical parameters d, and the individual likelihoods 



are combined by the methods described in Section 2.1 



2.3.1. The "Maximal Ignorance" Estimate 

The marginal posterior densities Ti2345(Tj, Cj), including the assumption of intrinsic alignment 
and the isophotal twist constraint (cf. equation [H]), are shown in Figure |l] for all models. This 
should be considered the "maximal ignorance" estimate of intrinsic shape, since as few restrictions 
as possible are placed on the dynamics. In each panel, oblate spheroids are at the right margin 
and prolate spheroids at the left. We rather conservatively allow values of c down to 0.35, but 
since no models with c > 0.8 are consistent, at even the 4<r level, with the observed e > 0.26, the 
upper part of the parameter space is truncated. White contours enclosing 68% and 95% of the 
total probability are the Bayesian equivalent of la and 2a error bars. 

The features to note in Figure p] are: (1) the strong similarity of the distributions for the inner 
4 annuli, with prominent peaks at T = 0; (2) the overall trend toward smaller c with increasing 
radius, consistent with the gradual radial increase in e; and (3) the significantly broader distribution 
in the outermost annulus, indicating an increase in T. The diagonal ridges in the distribution and 
the scalloped contours toward smaller c are numerical artifacts arising from undersampling of the 



sharply peaked densities -Pi(T, Cj, fij). Changing the outermost velocity on PA 225 (Sec. 2.2.3 ) 
causes a minuscule shift of the distributions for the inner four annuli toward smaller T and larger 
c, and in the fifth annuls flattens the secondary maximum near T = 0.2, c = 0.45 and truncates 
the weak appendage extending toward T = 1, c = 0.7. 

Figure |I| suggests that NGC 1700 is likely to be close to axisymmetric, at least in the inner 35" 
or so. Considering its central cuspiness, it is important to derive a quantitative constraint on the 
triaxiality T. By integrating the distributions in Figure [l] over c, we obtain the marginal densities 
for T shown in the first two panels of Figure [2|. Figure §a shows the differential distributions; the 
corresponding cumulative distributions are plotted in Figure The distributions for the inner 4 
annuli are so similar that one is naturally led to suspect that T may be constant over this range 
of radii. Indeed, plotting the integrated distributions tends somewhat to mask a strong preference 
in the joint distribution for constant T. This is demonstrated in Figure [2|c, in which we plot the 
differential probabilities for the change in triaxiality between annuli 1 and 2, 2 and 3, and 3 and 
4. All are sharply peaked at AT = 0. The strong similarity of these distributions and the small 
scatter in mean values are consistent with AT being drawn from a parent distribution that is a 
delta function at AT = 0, indicating constant T in the inner 4 annuli. 

If T is assumed to be the same in annuli 1-4, then this assumption can be included as "prior 



4 From here on we omit the L subscript on c, but it is still the axis ratio of the luminosity distribution. 
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Fig. 1. — "Maximal ignorance" estimate of the intrinsic shape of NGC 1700 in five radial zones. Each panel shows 
the marginal posterior density for the triaxiality Tj and short-to-long axis ratio Ci of the ith annulus, assuming that 
there is no intrinsic twist of the principal axes. Oblate spheroids are at the right margin and prolate spheroids are at 
the left. All velocity field models are used in this estimate, placing as few restrictions as possible on the dynamics. 
White contours enclose the 68% and 95% higest posterior density (HPD) regions. Note the strong similarity of the 
distributions for the inner four zones, the gradual tendency for greater flattening with increasing radius, and the 
marked change in the outermost annulus. 

information" (essentially by inserting factors of the form 5(Ti — Tj + i)J(Tj + i — Tj+2) • • • in equation 
p0|| ). The resulting cumulative distribution for the single value of T is shown in figure ^L. Upper 
limits on T at the la and 2a confidence levels are drawn as dotted lines. We obtain 

T < 0.16 for 2.2/T 1 kpc < r < 6.7/i _1 kpc (95% confidence). (4) 

A more intuitive but less precise statement is that the middle-to- long axis ratio b/a is greater than 
0.9 at > 95% confidence. 
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Fig. 2. — (a) Differential and (b) cumulative triaxiality distributions obtained by integrating the posterior densities 
in Figure |l| over c. The nearly identical distributions in the inner 4 annuli imply constant T over this region, as 
do (c) the differential probabilities for the change in triaxiality between adjacent annuli. Including the assumption 
of constant T in the analysis gives the cumulative distribution for the single value of T shown in (d) . Dashed lines 
indicate la and 2a upper limits. 



2.3.2. Statistical Effects of Orientation 

Because P(T, c) is a marginal density integrated over orientation Q, its value at a given (T, c) is 
proportional to the amount of solid angle over which a model of that shape will be seen to reproduce 
the observations (averaged over the other parameters). Consequently, there is a complementarity 
between constraints on shape and orientation. Shapes able to fit the data in many orientations are 
statistically favored while those requiring special orientations are not; orientations that are required 
by many shapes earn high marginal probabilities while those useful to only a few score low. This 
leads to the somewhat counterintuitive result that constraining the models to be at the same 
orientation has the effect of narrowing the marginal distribution in orientation while broadening 
that in intrinsic shape. This effect can be rigorously demonstrated for analogous one-dimensional 
problems. For NGC 1700's highly symmetric velocity field, triaxial models requiring lines of sight 
in or near the xz plane to suppress apparent minor axis rotation are made relatively more likely 
by the alignment constraint than axisymmetric models which can be more freely oriented. The 
resulting marginal densities -P(f2) thus show a strong preference for such orientations. 
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One should remember, however, that this estimate is made in isolation, nearly independent 
of knowledge of the rest of the elliptical galaxy population. That the data may imply a special 
orientation for a single object is hardly alarming. It is only in conjunction with other objects 
that the statistics are affected by the requirement that the observed sample must be consistent 
with a random drawing from an isotropic parent distribution. Since most observed ellipticals are 
known not to have strongly triaxial-like kinematics, this added constraint will push down the peak 
in the orientation distribution, eliminating triaxial models that require that special orientation 
and narrowing the marginal distribution for T. Methods for including such constraints have been 
developed (Bak & Statler 1997) and will be discussed fully and applied in a forthcoming paper (Bak 
& Statler 1998, in preparation). For now it suffices to acknowledge that including such constraints 
in the modeling of NGC 1700 would tighten the limits on T, so that the limit in equation (Q) can 
be regarded as conservative. 



2.3.3. The Dynamical Prior 

The maximal ignorance estimate intentionally places as few restrictions as possible on the 
orbit populations, because we have yet to understand how Nature chooses to populate orbits. Such 
knowledge may be forthcoming, however, either from observations or numerical simulations. Prior 
restrictions on the dynamics will alter the shape estimate. 

The infrequent occurrence of minor axis rotation in ellipticals seems to indicate that the "in- 
trinsic misalignment" between the angular momentum vector and the short axis tends to be small 
(Franx, Illingworth, & de Zeeuw 1991). It may be the case that net streaming motions are present 
only in the short-axis tube orbits. If all models with long-axis tube streaming are excluded from 
the fits, the result is the set of posterior densities shown in Figure ||[ The galaxy comes out more 
triaxial, since, with long-axis tube streaming disallowed, triaxial systems are more likely to be seen 
in orientations where the apparent minor axis rotation is zero. In Figure |]a we show the marginal 
densities for T alone. Those for the inner four annuli are not quite the same, but the differences, 
compared with Figure |2|a, are consistent with the larger amount of numerical noise in the restricted 
set of models. Adding the prior constraint that T is constant over this range of radii gives the 
cumulative distribution in Figure and a 95%-confidence upper limit of T < 0.22. Note that this 
is a change of only 0.06 from the maximal ignorance limit. 

By far the most important dynamical parameter is not the relative loading of long- and short- 
axis tubes, but the form of the v*(t) function, which gives the latitude dependence of the mean 
streaming velocity across the xz plane. Restricting the models to v*(t) types 1 and 4 or to types 2 
and 3 gives the marginal T probabilities shown in figure ||. (No restriction is placed on long-axis 
tube streaming.) The former case shows a marked preference for triaxial shapes over axisymmetric 
ones, the latter just the opposite. 

The reason for this difference is the way in which the models fit the position-angle dependence 
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Fig. 3. — Intrinsic shape estimate as in figure |l| but excluding all models with net streaming in long-axis tubes. 
The galaxy is inferred more likely to be triaxial, since triaxial systems forced to rotate only about their short axes 
are more likely to be seen in orientations where the velocity gradient along the apparent minor axis is zero. 

of the observed velocity field. SSC find that velocities on the sampled PAs are consistent with an 
almost exactly sinusoidal dependence at a given radius. To understand how this affects the models, 
we note first that the results for types 1 and 4 separately are nearly the same, as are those for 
types 2 and 3; this similarity indicates that the model velocity fields are dominated by short-axis 
tubes. In triaxial systems, the sequence of short-axis tubes terminates in the xz plane at a latitude 
cos" 1 T 1 / 2 up from the xy plane. The streamlines diverge to either side of this constriction (see fig. 
2b of S94a). The mean flow in the vicinity of the short axis will be seriously diluted, either by this 
divergence alone, or by additional geometrical factors if the line of sight is not close to the xz plane 
- unless either the flow velocity through the constriction is sufficiently large near the terminating 
latitude, or T is sufficiently small that the terminating latitude is itself close to the z axis. The 
type 2 and 3 models, where the flow velocity through the constriction falls off smoothly to zero 
with increasing latitude, rotate too slowly on the diagonal position angles relative to the major axis 
in projection, unless they are so close to oblate that the short-axis tubes extend nearly all the way 
to the z axis. Conversely, the type 1 and 4 models, where the velocity is assumed constant right 
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Fig. 4. — (a) Differential triaxiality distributions as in Figure fyi, but for the case with net rotation in long-axis 
tubes omitted, (b) Cumulative distribution for T assuming the value is constant in the inner 4 annuli. Dashed lines 
indicate la and 2<r upper limits. 
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Fig. 5. — Differential triaxiality distributions resulting from different choices of the form of the v*(t) function 
(Appendix B). (a) Types 1 and 4; (b) types 2 and 3. (See Section 2.3.3.) 



up to the terminating latitude, rotate too fast on the diagonal PAs unless the terminating latitude 
is rather low, meaning T is significantly above zero. 

Knowing the actual rate of occurrence of triaxiality in ellipticals has major consequences for the 
abundance of steep central cusps and/or massive black holes, and is an important test of the idea 
that central mass concentrations can destroy triaxiality. The above result implies that triaxiality 
constraints could be considerably strengthened if we could restrict, on physical grounds, the form of 
v*(t). It is, therefore, extremely important to obtain the v*(t) function from simulated ellipticals, 
e.g., from iV-body merger experiments (Dutta, Statler, &; Weil 1998). Furthermore, understanding 
how this function is determined physically by different mechanisms could make it possible, with a 
sufficiently large sample, to use present-day kinematics to constrain galaxy formation physics. 
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3. Mass-to-Light Ratio and Anisotropy 

Given the upper limit on T derived in the previous section, we can safely use axisymmetric 
models to estimate the mass-to-light ratio in NGC 1700. A wealth of techniques have been developed 
over the past several years to fit axisymmetric models to observations, of which two suffice for our 
purposes. 



3.1. Moment Models 

3.1.1. Method 

For axisymmetric systems in which the phase space distribution function / depends only on 
energy and the z component of angular momentum, the solution of the Jeans equations for the low 
order moments of / can be broken down into one-dimensional integrals along lines parallel to the 
symmetry axis and derivatives of these integrals (Satoh 1980; Binney, Davies, & Illingworth 1990, 
hereafter BDI). Efficient methods for constructing such models were developed by BDI and refined 
by van der Marel, Binney, & Davies (1990; hereafter MBD) and Magorrian & Binney (1994). Our 
implementation essentially parallels that of BDI and MBD, with minor modifications. 

The i?-band surface photometry is first interpolated onto an elliptic grid, spaced uniformly 
in position angle ip and logarithmically in elliptic radius m p = i?(cos 2 if) + sin 2 t/j/c 2 ) 1//2 ; here R is 
projected radius on the sky and we take c p = 0.7. Since we are assuming axisymmetry, PA twists 
are ignored. The photometry is extrapolated to small and large radii using separate r 1//4 -law fits 
to the inner and outer parts of the radial profile. 

The deprojection for a given inclination i starts with a parametric least-squares fit of a lumi- 
nosity density of the form 

PL{m) = — ^— j=. (5) 

In equation (|5|), a is a scale length and m is the spheroidal radius r(cos 2 # + sin 2 #/c 2 ) 1//2 , where 
c = [(c 2 — cos 2 £)/ sin i} 1 / 2 is the intrinsic axis ratio that projects to the apparent ratio c p at 
inclination i. This model projects to a surface brightness 

I(m p ) = tan -1 — , (6) 
m p m p 

which is a better fit to the NGC 1700 data than the Hernquist and Jaffe models used by MBD. The 
model is evaluated on a spheroidal grid, and the fit is then improved by iterating via Lucy's Method, 
giving a final luminosity model that is still axisymmetric but not spheroidal. Three iterations are 
sufficient to bring the photometric residuals down below 0.05 magnitudes. 

The mass density is found from 



p(r,9) = Tp L (r,e), 



(7) 
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Fig. 6. — The observed root-mean-square velocity profiles (points with error bars) are compared with predictions 
of the the best 2- integral constant-M/L model (dotted line) and a 2- integral model with a dark halo (dashed line). 
The constant-M/L model fails miserably while the halo model gives an acceptable fit. Both models are inclined at 
i — 90° (edge-on); models at lower inclinations are unable to match the major to minor axis dispersion ratio. 

where we can take the mass-to-light ratio T to be either constant or a function of m. We obtain the 
gravitational potential $ by a spherical harmonic expansion, calculate the mean-square velocities 
using BDI's equations (7) - (8), and project them onto the line of sight. We choose not to separate 
the ordered rotation from the random tangential velocities since this introduces more free param- 
eters into the models. Instead, we add the observed dispersions and mean velocities in quadrature 
and compare the resulting RMS with the models. 
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3.1.2. Results 

The dotted lines in figure |6| show the RMS velocity profiles for the best constant-M/L model. 
This model has i = 90° (edge-on) and has been fitted by eye to the data in the inner 10". The 
model clearly fails to reproduce the high RMS velocities observed beyond r e . Note also that for 
R < 10", the velocity on the major axis is slightly overestimated and that on the minor axis slightly 
underestimated. This mismatch only gets worse for flatter models at lower inclinations, since, to 
support the flattening, the tangential dispersion must increase at the expense of the radial and 
vertical dispersions (van der Marel 1991). We conclude that 2-integral models with constant M/L 
cannot fit the data. 

Two-integral models with variable M/L fare better, though it is difficult to constrain the 
details of the mass profile. To estimate the magnitude of the required M/L gradient, we adopt a 
simple law of the form 

T(m) = T [l + (m/r h ) b ] 1/b , (8) 

in which T is stratified approximately on density surfaces. For large b, this describes an M/L 
ratio that is constant at small radii, turns up abruptly at r^, and then increases linearly with 
radius. Since the surface brightness is still falling roughly as R~ 2 where the photometry ends, 
T ex r corresponds to what one expects from a standard isothermal dark halo. Fixing the value 
of b = 8 and minimizing x 2 with respect to To and Th results in the model shown as the dashed 
lines in figure § the fit is statistically adequate, with a reduced \ 2 of 1-16 (49 degrees of freedom). 
This model is inclined at 90°, since, as before, lower inclinations are unable to match the major to 
minor axis dispersion ratio. Making the halo spherical rather than spheroidal only exacerbates the 
problem, since then the galaxy cannot rely on the flattening of the potential to keep the luminosity 
density flat, and must become even more tangentially anisotropic. 

The best fit local and cumulative M/L profiles are shown in figure 0. We have converted to 
blue luminosity using a mean B — R of 1.7 (FIH). The central value of 5.5h agrees very well with 
Bender, Burstein, & Faber's (1992) result of 2.8 for h = 0.5, and the factor of ~ 2 cumulative 
increase out to a scale of ~ 10ft, -1 kpc is consistent with standard halo models. Thus 2-integral 
models with dark matter appear able to fit the data. 

This result, however, is less than entirely satisfactory for at least three reasons. First, the 
velocity field fitting of Section |2| implies flatter shapes and lower inclinations than obtained here. 
Second, the mere existence of a solution to the Jeans equations does not guarantee the existence 
of a non-negative distribution function /. And third, since mass can be traded off with anisotropy, 
3- integral models are essential to determine whether an M/L gradient is actually required by the 
data. 



3.2. Quadratic Programming Models 
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Fig. 7. — Mass-tolight ratios, in terms of local density ratios and enclosed quantities, plotted as functions of radius 
for the best 2-integral halo model shown in Figure ^[ 



Three-integral models constructed by quadratic programming (QP) score over moment models 
in that, since one solves explicitly for a positive definite phase space distribution function, there is 
no risk that an apparently successful model may require an unphysical negative density. 

The numerical approach used here is essentially the one described by Dejonghe et al. (1996). 
In brief, the procedure is as follows: first, for an assumed inclination i, the surface photometry is 
deprojected as described above using Lucy's method. The total mass distribution is assumed to 
have the same shape as the luminosity density, though not necessarily the same radial profile. This 
is accomplished by setting, for each term in the spherical harmonic expansion of the density, 



the case of constant M/L then corresponds to B = 0. The potential is computed from p by 
harmonic expansion. A Stackel fit is made to the potential, using the de Zeeuw & Lynden-Bell 
(1985) method as implemented by Dejonghe & de Zeeuw (1988). This fit is made only to define an 
approximate third integral ^3; wherever the potential is needed during the modeling, the original 
potential is used. The distribution function / is built up as a sum of components of Fricke type: 



3.2.1. Method 




(9) 




(10) 



p,g,n 



where n is an integer but p and q may be real, and 
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For components of this type, the velocity moments can be expressed in terms of algebraic functions 
of the coordinates and the potential. Models are constructed from a library of components with 
various combinations of p, q and n, generated by the direct product of p = 2 ... 10, q = . . . 5 and 
n = ... 4. The QP algorithm is used to compute the coefficients A^ qn that optimize the fit to 
the photometric data and the mean velocity and dispersion profiles. This is done by minimizing a 
X 2 -like variable of the form 

2 V - ^ [ Mobs,i P>mod,i 

x — 2_i Wi 1 

a \ A'modji 

where the //'s represent the photometric and kinematic data collectively and the sum is over all data 
points. While complete two-dimensional surface photometry is available, we choose for convenience 
to fit the photometry only at the points where there are kinematic measurements. This allows 
some of the models to come out a bit peanut-shaped, but the effect on the inferred mass profile is 
minimal. The weights Wi are in principle arbitrary, but are here assigned the value 1. 




3.2.2. Results 

The best constant- M/L 3- integral model, at an inclination of 90°, is shown in figure ||. The 
global value of M/Lb = 2.8h has been chosen to optimize the overall fit, i.e., without assigning 
larger weight to any particular region, though the higher S/N at smaller radii makes the model 
appear pinned to the central dispersion. The shallow slope of the RMS velocity profile cannot be 
reproduced; thus, allowing for anisotropy does not eliminate the need for an M/L gradient. In fact, 
the fit is even slightly worse than for the 2-integral moment model of § 3.1.2| . We ascribe this to 



the positivity requirement on / in the QP models, which is not applied in the moment models. If 
the positivity constraint is lifted, the same QP components can produce a fit with a \ 2 a factor 
of 3 lower. We take this as an indication that the constant-M/L moment model actually does not 
correpsond to a nonnegative distribution function. 

Allowing for an M/L gradient improves the fit substantially. The dashed lines in figure || show 
a model with i = 90°, p = 1, B = \.Tlh~ x kpc -1 , and To = 5.1h (B band, in solar units), compared 
with the observed rotation (squares) and dispersion (crosses) profiles. The fit to the kinematics 
alone, ignoring the photometry, is passable, with x 2 l v = 2.03 for v = 108 degrees of freedom. 
The discrepancy at large radii on PA 315 is not especially worrisome since this is where isophotal 
and kinematic twisting sets in, and axisymmetric models should not be expected to reproduce the 
detailed velocity field. 

The local and cumulative M/L profiles for this model are shown as the dashed lines in figure 
10. The central values agree very well with those shown in Figure ^ for the 2-integral model, and 



appear to be insensitive to assumptions about the functional form (i.e., equation |8| vs. equation 
|J.) At large r the local M/L gradient is roughly a factor 1.5 smaller than that obtained from the 
2-integral model. This translates to a ~ 20% difference in the cumulative M/L gradient because 
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Fig. 8. — The best three-integral constant- M/L model compared with the observed root-mean-square velocity 
profiles, as in figure ^| The added freedom in the distribution is insufficient to reproduce the shallow profile at large 
radii. 



of the steepness of the mass profile. 

Obtaining a strict lower limit on the M/L gradient, while desirable, is not especially practical. 
A shallower gradient can be obtained at the expense of a more pathological distribution function 
at low binding energy, corresponding to distortions of the LOSVD at radii for which there are no 
data. As there is no objective criterion for accepting or rejecting such models, we have rather 
subjectively chosen the parameters B and p so as not to overly tax the QP algorithm; thus the 
distribution function is not on the verge of going negative, and the anisotropy is not extreme. 



The latter property is illustrated in the top two panels of figure |11| which show the ratios of the 
unprojected dispersions ctr, a^, and a z in the meridional plane. The velocity ellipsoids are generally 
elongated in the tangential direction, and (j z /or is within 10% of unity nearly everywhere. The 
model becomes more tangential at larger R; but remember that the kinematic data end near 70", 
so the velocity distribution beyond this point is not well constrained. 
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Fig. 9. — Three-integral, variable-M/L models that adequately fit the kinematic data, compared with the observed 
mean velocities (squares) and dispersions (crosses). Models have inclinations of 90° (dashed lines) and 60° (dotted 
lines). The 60° model has been scaled up by 3% from the formal best fit to better reproduce the central dispersion. 



Though we have not checked explicitly, it seems virtually certain that for any inclination % > 50° 
(for which c > 0.35) it is possible to find many values of B and p that produce an acceptable fit. For 
the shape fits of Section |2|, the posterior densities in figures |l| and |3] give a most-probable intrinsic 
flattening c between 0.55 and 0.65, which would correspond, for an axisymmetric system, to i « 60°. 
A QP model at this inclination, using the same B and p as the i = 90° model above, is shown as 
the dotted lines in figure ||. The formal best fit has been scaled up by 3% in velocity, or 6% in mass, 
to improve the fit to the central dispersion. The local and cumulative M/L profiles are plotted as 



the dotted lines in figure 10, and the ratios of the principal dispersions are shown in the bottom 



two panels in figure 11. This rather flat model is tangential at all radii, and at large R (past where 
the data end) becomes radially quite cold. Such a configuration is not at all impossible, especially 
if, as SSC suggest, the main body of the galaxy is surrounded by a differentially precessing disk 
acquired through a merger or other accretion event. 
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Fig. 10. — Local and cumulative mass-to-light ratios, as in Figure [?], for the i — 90° (dashed lines) model and the 
i = 60° model (dotted lines) shown in Figure 0. 



4. Discussion 

By modeling its velocity field and surface photometry, we have found that NGC 1700 is nearly 
axisymmetric and oblate. We can say with better than 95% confidence that if it were viewed down 
its short axis, it would appear as at most an El at radii <^ 36" (6.7/i _1 kpc). Beyond this point, if 
the density surfaces remain aligned, then photometric and kinematic twists indicate that it becomes 
significantly flatter and more triaxial. 

The dominant source of systematic error in the determination of intrinsic shape is our lack 
of prior knowledge about the internal dynamical structure, specifically, how we should expect the 
mean streaming velocity to vary away from the symmetry planes. Axisymmetry is strongly favored 
if the velocity declines away from the xy plane; triaxiality is preferred if it doesn't. It is therefore 
important to understand how the latitude dependence of the mean streaming may be determined 
by physics. Preliminary results of an analysis of group merger simulations (Dutta, Statler, & Weil 
1998), though limited by small-number statistics, suggest that the mean velocity does decline away 
from the xy plane in objects produced by this mechanism. It remains to be seen whether different 
formation scenarios would produce different behavior. 

NGC 1700 falls into the family of steep-cusped ellipticals delineated by WFPC photometry, 
standing out, in fact, as its most luminous member (so far). In general, the steep-cusped galaxies 
show characteristics one would expect of axisymmetric systems, such as disky isophotes and rel- 
atively rapid rotation. The near axisymmetry we have found is consistent with Merritt's (1997) 
result that orbital stochasticity prevents the existence of self-consistent triaxial equilibria for the 
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Fig. 11. — Ratios of the internal (i.e., not projected) velocity dispersions for the (top) i = 90° and (bottom) i = 60° 
models of Figure []. The ratios of the (left) azimuthal and (right) vertical components of dispersion to the radial 
component are contoured in the meridional plane. Kinematic data end at around 70". 



steep (r -2 ) cusped Jaffe profile unless T 0.4, T ^ 0.9, or 0.8. It is tempting to interpret the 
change of shape beyond 40" as marking the limit out to which triaxiality has been eaten away by 
stochastic scattering of box orbits. But it is far from clear that the cusp in NGC 1700 is the cause 
of this structure. As discussed by SSC, there is reason to believe that the observed photometric and 
kinematic twists arise from an outer accreted stellar disk or ring which is incompletely phase-mixed 
and therefore intrinsically misaligned with the rest of the galaxy. In this case the dynamical time at 
the location of the twist, times a factor of order a few, indicates the time since the last significant 
merging event. The dynamical age of 2 to 4/i _1 Gyr (SSC) is corroborated by Whitmore et al. 
(1997), who have observed NGC 1700's globular clusters and find colors and magnitudes consistent 
with a comparably old cluster population. 

We have also found, using axisymmetric 2-integral moment models and 3-integral quadratic 
programming models, that NGC 1700 must have a radially increasing mass to light ratio. No 
constant-M/L model can fit the data. The cumulative (enclosed) M/Lb increases from a central 
value of ~ 5h to ~ Wh (in solar units) at r = 12/i _1 kpc. The case for dark matter in NGC 
1700 is strong; even allowing for varying anisotropy, the best constant-M/L model fits the RMS 
velocity profiles at no better than x 2 l v = 5.74 (for v = 50). In this sense the evidence of dark 
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matter is comparable to that in NGC 2434 (Rix et al. 1997); in fact, NGC 1700 is arguably a 
more compelling case since only spherical models have been applied to NGC 2434. Other systems 
suspected on stellar dynamical grounds of harboring substantial dark halos, such as NGC 1399, 
4374, 4472, 5813, 7626, 7796, and IC 4296 (Saglia, Bertin, & Stiavelli 1992, Saglia et al. 1993, 
Bertin et al. 1994) are statistically less convincing. NGC 1700 appears to represent the strongest 
stellar dynamical evidence to date for dark halos in ellipticals. 

Bertin et al. (1994) point out that, of the ellipticals listed above with stellar dynamical indica- 
tions of dark matter, most turn out to be known bright X-ray sources. Unfortunately, NGC 1700 
seems never to have been observed by an X-ray telescope, though ROSAT programs are pending. 
NGC 1700 bears a more-than-passing resemblance to NGC 7626 (Saglia, Bertin, & Stiavelli 1992), 
having a very similar luminosity, central dispersion, effective radius, and distance. However, NGC 
7626 is 2 orders of magnitude more powerful in the radio (Bender et al. 1989), and also resides 
in a rich cluster while NGC 1700 has only one nearby companion and lies among small groups. 
It will be very interesting to see whether these two apparently similar ellipticals, in very different 
environments, turn out to be X-ray twins. 

TSS acknowledges support of this work from NASA Astrophysical Theory Program Grant 
NAG5-3050 and NSF CAREER grant AST-9703036. We thank the director and staff of the Multiple 
Mirror Telescope Observatory for their generous assistance and allocations of time to this project, 
and Trevor Ponman and Joe Shields for timely information. 



Appendix A: Combining Shape Estimates for Multiple Annuli 

We begin by introducing the notation Pijk...(Ti,Ci, Qj) for the marginal probability density for 
the shape and orientation of annulus i, making use of the data from annuli i,j,k,.... Suppose 
that we have obtained Pi(T\, ci, Ox) and P<i(Ti^ C2> ^2) for annuli 1 and 2 individually (c/. equation 
0). (Here both T and c refer to the light distribution; we omit the L subscripts for clarity.) We 
assume that the luminosity distribution is intrinsically aligned (Q,\ = O2 = O) and require that the 
observed position angle difference 812 ±c<5 12 be reproduced. Let Lg 12 (Ti, T2, f2) be the likelihood of 
the observed twist given the pair of intrinsic shapes (Ti, ci) and (T2, C2), and orientation Q: 

L 5l2 (T 1 ,T 2 ,n) = expj j. (13) 

The axis ratios c\ and c 2 do not appear in either the likelihood L$ 12 or the predicted twist 5f2° d 
because in the quasi-local approximation the major axis position angle in annulus i, 

PAi = I tan" 1 —— 2 T ^ 6 ^ 2 ^ (14) 

2 |7i(sin 4>e — cos 2 4>e cos 2 6e) — sin 4>e\ 
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is independent of c;. The joint probability that the two shapes fit the local photometry and 
kinematics and the PA twist is given by 

P(T 1 ,c 1 ,T 2 ,c 2 ) = [ dnP 1 (T 1 ,c 1 ,n)P 2 (T 2 ,c 2 ,n)L Sl2 (T 1 ,T 2 ,n), (15) 

but this is a somewhat unwieldy quantity. Of more use is the marginal probability for, say, annulus 
1 alone: 

iWi.ci) = / dT 2 dc 2 P(T 1 ,c 1 ,T 2 ,c 2 ) = J dOP^Tx, 01,0)^12(^,0), (16) 
where the "twist kernel" K\ 2 is given by 

K 12 (T 1 ,Q) = f dT 2 L 5l2 {T u T 2 ,^)P 2 {T 2 ,9), (17) 
and P 2 (T 2 ,Q) is the marginal distribution integrated over c 2 , i.e., 

p 2 (t 2 ,q) = J dc 2 p 2 (T 2 , C2 ,n). (is) 

The posterior -Pi2(7i,ci) is the shape estimate for annulus 1 considering all of the data in both 
annuli, but regardless of the actual shape of annulus 2. 

In our case the Pi(Ti,a,£li) are obtained from the dynamical models, so the Tj are mass 
triaxialities. If one is unwilling to assume a relation between the triaxialities of mass and luminosity, 
one can still combine the annular fits, assuming only that the principal axes are intrinsically aligned 
and foregoing the isophotal twist constraint, by letting a^ 12 —> oo, which gives L$ 12 = 1 by equation 

©• 

The generalization to N annuli uses the measured PA gradients 6± 2 , S 2 3, . . . , Sn—i,n, as follows: 
Pi2...N(Ti,Ci,n) = PiiTi^tyKflL^Ti, n)Kf^ hmtN (T if il), (19) 
where the interior and exterior kernels are given by 

#§-i,...,i (r*,0) = J dT i ^P i ^(T i . 1 ,Q)L 5l ^ l (T l -i,T l ,n) 

J dTi^Pi^Ti^tyLs^.^iTi^T^un)--- (20a) 

J dT 1 P 1 (T 1 ,Q)L Sl . 2 (T 1 ,T 2 ,n); 

K?? +l _ N (T u n) = JdT i+1 P i+1 (T i+1 ,n)L s . i+1 (T i ,T i+1 ,n) 

J dT i+2 P i+2 (T i+2 ,n)L Si+1 . +2 (T i+1 ,T i+2 ,n)--- (20b) 

J dT N P N (T N , Q)Ls N _ 1JV (T/v_i, T/v, O). 

As before, the shape estimate for annulus i, using all of the data but regardless of the shape at 
other radii, is the marginal distribution Pi 2 ...N(Ti,Ci) = J dfl P\ 2 ,,,]y(Ti,Ci,Q). The construction 
of equations (|l9|) and (|20|) guarantees that the marginal distribution in orientation, -Pi2...7v(0) = 
/ dTi dciPi 2 ,..N(Ti,Ci,£l), is independent of i, as it must be given the assumption of intrinsic align- 
ment. 
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Appendix B: The Contrast and v*(t) Functions 

In the models used in Section H properties of the phase space distribution function are sub- 
sumed into the "similar flow" ansatz, a scalar constant C, and a function of one variable v*(t). The 
last two describe the mean velocity across the xz plane on one fiducial spherical or ellipsoidal shell, 
which in turn determines the velocity field over the whole shell once T and the luminosity density 
are specified. 

The "contrast" C is defined as the ratio of the y component of the mean velocity on the x axis 
to that on the z axis, on the fiducial shell. One may think of this as the ratio of the maximum 
mean streaming velocities produced by short- axis and long- axis tubes. Models lacking short-axis 
or long-axis tube streaming have C = or C = oo, respectively. Models are run with values of C 
both independent of intrinsic shape and given by the four functions of (T, c) shown in Figure 2 of 
S94c. 

The function v*(t) gives the angular dependence of the mean velocity across the xz plane on 
the fiducial shell. The variable t is a rescaled polar angle, defined so that t = corresponds to the 
x axis, t = 2 to the z axis, and t = 1 to the locus dividing the short-axis from the long-axis tubes. 
For spherical shells, 

J 2 -^, ^sin- 1 ^, 

for ellipsoidal shells the relation between t and 9 is more complicated, and we refer the reader to 
Section 3.1 of S94c. By definition u*(0) = C and v*(2) = 1. In the intervals [0, 1) and (1, 2], v*(t)is 
assumed to be either piecewise-constant or piecewise linear, in the latter case falling to zero at 
t = 1. This gives four possible forms for v*(t): Type 1, constant in both intervals; Type 2, linear 
in both intervals; Type 3, linear in [0, 1) and constant in (1,2]; and Type 4, constant in [0, 1) and 
linear in (1,2]. 
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